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ABSTRACT 

I present empirical measurements of the rate of relaxation in A-body simulations 
of stable spherical systems and distinguish two separate types of relaxation: energy 
diffusion that is largely independent of particle mass, and energy exchange between 
particles of differing masses. While diffusion is generally regarded as a Fokker-Planck 
process, it can equivalently be viewed as the consequence of collective oscillations that 
are driven by shot noise. Empirical diffusion rates scale as N~^ in inhomogeneous 
models, in agreement with Fokker-Planck predictions, but collective effects cause re¬ 
laxation to scale more nearly as in the special case of a uniform sphere. I use 

four different methods to compute the gravitational field, and a 100-fold range in the 
numbers of particles in each case. I find the rate at which energy is exchanged between 
particles of differing masses does not depend at all on the force determination method, 
but I do find the energy diffusion rate is marginally lower when a field method is used. 
The relaxation rate in 3D is virtually independent of the method used because it is 
dominated by distant encounters; any method to estimate the gravitational field that 
correctly captures the contributions from distant particles must also capture their 
statistical fluctuations and the collective modes they drive. 

Key words: Galaxies: kinematics and dynamics — methods: numerical 


1 INTRODUCTION 

The topic of relaxation driven by stellar encounters in 
star systems has a l ong and distinguished his tory {e.g. 
IChandrasekhaj Il94ll : iBinnev fc Trern^aind l2008l l and has 
important im plications for the evolution o f star cluster s 
llSpitzerj|l987ll and of active galactic nuclei (lMerritt|l2013ll . 
However, it is widely believed that relaxation through star- 
star encounters occurs at a negligible rate in the bulk 
of galaxies, and therefore A-body simulations of galaxies 
should mimic this collisionless property. The study presented 
here focuses on just one small aspect of the general problem 
of relaxation in simulations. 

Although the computational power available to simula¬ 
tors has risen steadily over time, calculations with billions of 
particles are still not routinely possible and most simulations 
employ fewer, more massive particles. A rough estimate of 
the relaxation t ime treiax in a simulation of N equal mass 
particles is (e.o. iBinnev fc Tremainj|2008ll 

trelax/tcross O.lA/lnA, (1) 

with fcross being a typical crossing time. This, as other more 
precise expressions for Fokker-Planck diffusion, contains the 
Coulomb logarithm that arises from integration over impact 
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parameters, implying that every decade of impact parame¬ 
ter makes an equal contribution to the integral. It should be 
noted that the approximate fo rmula (1) applie s for pressure- 
supported systems, not discs llSellwoodll2013l l. and neglects 
collective effects. Furthermore, it was derived for nominally 
point mass particles, since the argument of the logarithm 
comes from the ratio of the system half-mass radius to the 
scale on which scattering causes large deflections, whose con¬ 
tributions would otherwise be overestimated. It therefore 
reflects spatial resolution, which is usually determined by 
other factors, such as particle softening or grid resolution 
in A-body codes. Therefore, treiax should be simply propor¬ 
tional to A when resolution is held fixed. 


Numerical methods used to compute the gravitational 
field in simulations that aspire to be collisionless fall into 
three broad categories. The most popular direct method is 
some type of tree code {e.g. iBarnes fc Hut|[l98^ : [s pringell 
I 2 OO 5 II that effectively sums the attraction of every par¬ 
ticle pair, with a softening kernel to limit the magni¬ 
tude of the acceleratio n at short range. The far more ef- 
ficient ( ISellwoodI 20l4l par ticle-mesh methods (PM, e.g. 
iHocknev fc EastwoodI Il98lll determine the gravitational 
field on a raster of points that has some appropriate ge¬ 
ometry; forces at the actual particle positions are com¬ 
puted by interpolation between grid points. Finally, the 
least popular are field methods {e.g. IClutton-BrockI [19721 : 
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iHernauist fc Qstrikeil[l9^ 1 that expand the density and po¬ 
tential in a basis set of functions that should be chosen such 
that truncating the expansion at low order yields an ade¬ 
quate approximation to the field. _ 

Wein berg, in a series of papers (IWeinberd 


I 999 I: Hollev- Bockelmann, Weinberg fc Katj 20051 : 


Weinberg fc Katzn2007al ^ has argued that field methods 


are inherently superior to other Wbody methods in their 
ability to hide the lumpiness o f the potential from a set 
of point masses. In particular, IWeinberg fc Kat3 l|2007bl 'l 
assert that field methods have “relaxation times ... orders of 
magnitude longer” than in tree codes. This claim was based 
on a lengthy calculation using Hamiltonian mechanics that 
I review below. 

Since a well-chosen basis can yield an adequate approx¬ 
imation to the total field from a small number of terms, it 
may seem reasonable to expect the resulting potential to be 
smoother than that computed by other methods. However, 
this argument is beguiling for the following reasons. Only the 
lowest order monopole term has a large value about which 
shot noise fluctuations are small, while the amplitudes of 
the aspherical terms, which are oscillatory, depend on the 
almost-complete cancellation of the N contributions, and 
are entirely noise-driven in a spherical model, for example. 
It is also true that the number of values that define the 
potential in PM codes is rigrid and, since each mesh point 
typically hosts ~ particles, each separate value will 

be subject to a greater degree of shot noise. In direct meth¬ 
ods, there are N separate contributions to the field. But 
note that the potential is the double integral of the density; 
in effect, the potential kernel, which is monotonic and has 
infinite range, implies the field at each point is the sum of 
contributions from all the sources. Thus potential variations 
are far smoother than could be supported by an arbitrary 
fun ction defined by the s a me nu mber of values. 

IHernauist fc Barnes! (Il990l . hereafter HB90) used a 
King model for an experimental comparison between the 
relaxation rates in simulations when the gravitational 
field was determined by three different methods, and 
IHernauist fc Qstrikeil (1 19921 ') extended their results to in¬ 
clude a field method. They found that the rate of energy 
diffusion of particles was only mildly affected by the method 
used. While their evidence was quite strong, they generally 
employed only 4096 particles and gravity softeni ng spoiled 
the eq uilibrium of some of their initial models. ISellwoodI 
1 I 2 OO 8 II showed that grid and field methods performed equally 
well in the specific problem of bar-halo friction. 

Here I present a more general study of fully self- 
consistent equilibrium spheres that uses two distinct mea¬ 
sures of the “relaxation” rate: the energy diffusion rate re¬ 
ported by HB90, and a measure of the rate of energy ex¬ 
change between particle species of differing masses. The for¬ 
mer measure includes all sources of relaxation, especially col¬ 
lective effects, while the latter is a more direct consequence 
of two-body encounters. 


All three are spheres with ergodic (isotropic) distribution 
functions (DFs) that are therefore stable equilibria. 

The mass models are: 


(a) Hernquist model IHernauist! (ll99Cll) developed a simple 
spherical model having the centrally cusped density pro¬ 
file 


p(r) 


Ma 

27rr(r -|- a)® ’ 


( 2 ) 


which has the potential 4>(r) = —GM/{a + r). Here a 
is a length scale and M the finite total mass integrated 
to infinity. Hernquist also gave the equilibrium isotropic 
DF for this density and potential. I truncate this model 
so that no particle has sufficient energy to stray be¬ 
yond r = 10a, which causes the density profile to taper 
smoothly to zero at that radius, and discards ~ 26% 
of the mass, but leaves the density profile as given by 
eq. © over the r ange 0 ^ r/a < 5. 

(b) Plummer spfeere lPlummerl (Il91ll) introduced one of the 
most widely used spherical mass models in astronomy. 
It has the cored density profile 


p{r) 


3M 


(1 + ® ) 


2 \- 5/2 


( 3 ) 


where x = r/a and M is the total mass. The potential 
is <I>(r) = —(GM/a)(l -|- and the isotropic DF 

is that of a polytrope of index 5 (see BT08). Applying 
an energy truncation so that no particle passes outside 
r = 10a discards only ~ 3.4% of the mass in this model 
with its lower density envelope. 

(c) Uniform sphere wit h density pp = 3M/(47ra^)j where 
a is the outer radius. IPolvachenko fc ShukhmM (Il979f) 
give the isotropic DF for a homogeneous sphere with 
a sharp outer boundary. The harmonic potential in the 
interior of this unusual model implies that all particles 
have the same orbital frequencies; this model, therefore, 
affords a dramatic illustration of the role of collective 
effects. 


In all cases, I employ two equally numerous sets of parti¬ 
cles drawn from the DF: the masses of particles, mi = p.imf 
are such that those of one sample have pi = 9 and the 
other have pi = 1, and m* is chosen such that the com- 
bined density profile is t hat g iven by the above expressions. 
IPebattista fc SellwoodI (l2000l . appendix A) describe an op¬ 
timal method of drawing particle coordinates from a DF in 
such a way as to redu ce shot noise in the distribution of ener¬ 
gies. ISellwoo3 ll2014h reports results that show the material 
benefit of this strategy. 

The models are evolved for 100 dynamical times, where 
to = {a^, with a timestep of 0.02to for the uni¬ 
form sphere and Plummer models. The basic step for the 
Hernquist model is 0.0125 but time steps are increased, in 
this case only, by four successive factors of 2 at appropriate 
radii. I save the instantaneous energy of a representative set 
of particles (2 x 10'* or all, whichever is the less) after every 
dynamical time. 


2 MODELS 

In order to illustrate the importance of collective effects, I 
here report measurements in three different mass models. 


3 METHODS 

I employ four different numerical methods to determine the 
gravitational field from the particles: 
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(i) BHT A tree code th a t use s the scheme hrst pro¬ 
posed bv iBarnes &: HutI (ligsd l. I include dipole terms 
and particle groupings are opened when they subtend 
an angle > 6'niax = 0.5 radians. Forces are softened 
at short r a nge o nly using the kernel advocated by 
iMonaghlinl (Il992l ~l. with e = 0.05a for the uniform 
sphere and e = O.lo for the other two models. 

(ii) S3D A spherical grid, which is a hybrid PM expan¬ 
sion method that expands the non-spherically symmet- 
ric part in surface harmonics on a set of spherical shells 
(lMcGlvnnlll984 [Sellwoodll2003h . Force discontinuities, 
which arise when particle radii cross, are eliminated 
by adopting linear interpolation between radial shells. 
I tabulate the expansion coefficients at 100 logarith¬ 
mically spaced radii, and expand in azimuth up to 

^max — 8. 

(iii) SFP A held method that e mploys a biorthonor¬ 
mal set of basis functio ns l|Clutton-Brockl 1 19721 : 
iHernauist fc Ostrikeil Il99'3 l. which I name SFP (for 
smooth held-particle) but is also known as SCF. I ex¬ 
pand in azimuth up to Imax = 8 and employ radial 
functions up to Umax = 10. I use this method for the 
Plummer and Hernquist models only. 

(iv) C3D I do not use the SFP method for the uniform 

sphere, b ut employ a P M method that uses 3D Carte¬ 
sian grid Jjame and I set a = 50 mesh spaces. 

The grid has 129® points, except for the lowest N case 
where it was 257® points in order to allow plenty of 
room for expansion of the particle distribution. Linear 
interpolation results in the inter-particle force given in 
ISellwood fc Merritd (| 19941 . appendix), which is well ap- 
proxim ated by cubic d ensity kernel with e « 1.8 grid 
spaces (ISellwoodll2014h . 

More de tails of all thes e methods are given in the on-line 
manual (|Sellwoodll2014l L 

The choices of numerical parameters in each case are 
somewhat arbitrary, but values are typical of those used 
in practice and are not varied as the particle number is 
changed. The rate of relaxation will depend only weakly 
on spatial resolution, since only short range scattering is af¬ 
fected by changes to the softening length, for example, lead¬ 
ing to a slight change in the value of the Coulomb logarithm. 

Truncating the expansion at low order in field methods 
smooths the mass distribution, and more aggressive trunca¬ 
tion will give rise to a smoother potential, which must reduce 
the relaxation rate; e.g. in the extreme case of a single term, 
the particles will be moving in an almost fixed potential, 
and no evolution and little relaxation could occur. Since the 
purpose of simulations is to compute the self-consistent evo¬ 
lution as the density changes, I include sufficient terms to be 
able to follow changes that might be expected in an evolving 
model. 

The tree code uses explicit particle softening, while 
short-range forces are implicitly smoothed in the Cartesian 
grid. These methods therefore do not yield the exact New¬ 
tonian potential of the mass distribution. In order to ensure 
that the initial model is in equilibrium in the tree code, I 
use the largest N simulation of each type to tabulate the dif¬ 
ference between the spherically averaged central attraction 
of the particles at t = 0 and the analytic central attrac¬ 
tion at a ID array of points; I then interpolate from this 


table a supplementary central attraction that is added to 
the tree-determined force on each particle before its motion 
is advanced. I use a similar procedure for the 3D Cartesian 
grid, but in order to avoid shot noise in the (very small) 
permanent part of central attraction, the numerical force is 
computed from a smooth density distribution assigned to 
the grid. These generally small corrections are not needed 
for S3D or SFP methods. 

With these fixed correction terms in C3D and BHT, 
and with a coordinate centre for the S3D and SFP methods, 
it is important to ensure that the initial set of particles is at 
rest and centred in this coordinate frame. After creating the 
initial set of particles, I therefore adjust the positions and 
speeds by a small amount to ensure the centre of mass is 
at the coordinate centre and the model has no net momen¬ 
tum. An alternative strategy that would achieve the same 
outcome would be to shift the grid (or expansion) centre at 
frequent intervals. I also experimented with inserting mirror 
pairs of particles having coordinates (x, v) and x, —v), 
in a step towards a quiet start, but this strategy had the 
undesirable (for this study) effects of eliminating any lop¬ 
sided, and emphasizing the bi-symmetric, contributions to 
the total field. 

The instantaneous energy of the ith particle is approx¬ 
imately 

Si = miEi{t) = TUi [>l>(a:i) -I- \vf\ , (4) 

where Ei is the specific energy of the particle, or energy 
per unit mass, $ is the estimated gravitational potential at 
the particle position, Xi{t), and Vi{t) is the scalar speed 
of the particle. This definition is exact only for a particle of 
inhnitesimal mass, since a hnite mass particle contributes to 
"F. The energy required to disperse a system of gravitating 
particles, the total energy, is Stot = T + W, where T = 
'^i^rnivf and W = i miF) Xi), and the W term is 
halved because the summation over the 4> values includes 
every pair of particles twice. The total energy, dehned this 
way is very well conserved in the simulations, whereas the 
sum of the energies (eq. |T|, . Si = T + 2W, is clearly not 

the total energy, and is not conserved. 

In order to test for virial equilibrium of an Ai-body 
simulation, it is better to measure the virial of Clausius, 
Wc = '^^miXi ■ ai, since the inter-particle forces are not 
perfectly Newtonian, especially at short range. It is easy to 
show that Wc = W for precisely inverse-square law accel¬ 
erations. The particle distribution is in equilibrium when 
2T = |lTc|. 

It is convenient to chose units such that G = M = 
a = 1, so that the dynamical time to — {a ^= 1, 
for example. A convenient scaling to physical units for the 
inhomogeneous models is to choose a = 3 kpc and to = 
10 Myr, which implies M = 6 x M© and velocities scale 
as {GM ~ 293 km/s. 


4 RESULTS 

Figure [T] shows the final {t = lOOto) density proHles in all 
27 simulations. The three panels show the different mass 
models. Within each panel the curves drawn in red are from 
the lowest = 4 x 10®, those in green employ 4 x lO'^ 
particles, and those in blue 4 x 10®, with N/2 in each mass 
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Figure 1. The final density profiles (coloured lines) measured from the different particle species that, in all 27 simulations, contribute 
10% and 90% of the total density as indicated by the two groups of lines. In each panel, the colour indicates the number of particles: red 
is for TV = 4 X 10^, green is for TV = 4 X 10^ and blue is for TV = 4 X 10®. The different methods for computing the field are distinguished 
by the line styles: solid is for S3D, dotted is for BHT, and dashed is for either SFP (Hernquist and Plummer) or C3D (uniform sphere). 
The black dotted curves indicate the expected density profile of each mass species; no particles in the left and middle panels started with 
enough energy to pass beyond r = 10a. 


species so that the initial densities from each species differ 
by factors of nine, as shown. The field is determined on the 
S3D grid for the full-drawn curves, the BHT method for the 
dotted curves, while the dashed curves are from either the 
SFP, for the Hernquist and Plummer models, or C3D for 
the uniform sphere. 

The density profiles of the models evolve insignificantly 
for the largest TV (blue curves), confirming that the mod¬ 
els are stable equilibria. As expected, relaxation drives the 
greatest changes in the simulations with the smallest N (red 
curves) where, in a few cases, there are hints of some slight 
segregation of the particles of different masses over the time 
interval computed. 

4.1 Collective modes 

Figure [2] shows the time evolution of the virial ratio, T/|H4| 
in all 27 simulations reported here. The three panels show 
the different mass models and the colours and line styles are 
used to distinguish the particle number and field determi¬ 
nation method, as in Fig. [T] 

It is clear from Figure [2] that the virial ratio T/|Wc| re¬ 
mains close to i for the duration of all these simulations, and 
the fluctuations around this ratio are largest for the lowest 
TV (red curves) and smallest for the highest N (blue curves). 
The fluctuations are aperiodic for the inhomogeneous mod¬ 
els, but are periodic in the uniform sphere, where they have 
very nearly the same period, 27rto, in all nine simulations. 
Furthermore, the initial behaviour of the curves for models 
with the same TV (colour) is very similar for the different 
field determination methods, because the models were set¬ 
up using particles having the same positions and velocities. 
The red curves diverge quite noticeably, but the blue curves 
for the SCF and S3D methods remain barely distinguishable 
to the end; a slightly different value of H4 arises in the BHT 


code because the forces differ due to softening and corrective 
terms, which is the reason the dotted curves remain distinct, 
particularly in the Hernquist model. 

It is helpful to think of these fluctuations as the con¬ 
sequences of collective modes that are driven by shot-noise 
in the finite number of particles, as was recognized long ago 
bv iRostoker fc Rosenblu^ (Il96d ) in the context of collision- 
less plasmas, and has been studied in gravitating systems by 
IWeinbei^ (1 19981 ) and others. The simplifying concept here 
is to distinguish the ideal collisionless system, which has a 
set of neutral and/or damped modes, from the noise spec¬ 
trum arising from the particles that excites the modes. In 
principle, this approach would enable the evolution to be cal¬ 
culated by perturbation theory. Naturally, such an analysis 
should not differ in its predictions from the direct evolu¬ 
tion of an TV-body system composed of equal mass particles, 
which is exactly what the simulations compute. The results 
presented below reveal that relaxation is largely independent 
of particle mass, indicating that global modes are the domi¬ 
nant relaxation mechanism even when not every particle has 
the same mass. 

The amplitudes of the modes should scale with the shot 
noise, and the rms scatter in T/\Wc\ about the value 0.5 
indeed scales approximately as . The modes have es¬ 

sentially the same frequency in the uniform sphere where 
they appear to be almost undamped. Modes in collisionless 
systems can be damped only through resonant exchanges 
between the mode and the particles, and conditions for res¬ 
onant damping are highly unfavorable because all particles 
have the same frequencies in this harmonic potential. 

The aperiodic fluctuations in the inhomogeneous mod¬ 
els (left and middle panels) reflect the broader range of fre¬ 
quencies among the particles in these models. Collective os¬ 
cillations that are excited by shot noise in these cases are 
almost certainly quickly damped at resonances, yet the am- 
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Figure 2. The time evolution of the virial ratio in all 27 simulations. As for Fig. [T] the colour indicates the number of particles while 
the different methods for computing the field are distinguished by the line styles. 



Figure 3. Results from a run of a Hernquist model with N = 
2 X 10^ particles of each mass when forces are determined using 
the field method (SFP). The left panel shows the mean square 
change of energy of the particles, while the right panel shows 
the evolution of the mean specific energy (per unit mass) of the 
particles of each mass species. 


plitudes fluctuate greatly with neither a clear decaying, nor 
a growing, trend. 


4.2 Energy diffusion 

The potential in a hypothetical simulation with infinitely 
many particles would be smooth and steady, assuming the 
model is a stable equilibrium, and the specific energy E of 
each particle would be conserved. Therefore the time evolu¬ 
tion of the rms changes in E is one convenient measure of 
relaxation in a system with finite N. 

Figure [3] presents a typical set of energy measurements 
from the particles in a simulation; values for heavy particles 
are drawn in red, while green is used for the light particles. 
The left panel gives the value of {[Ei{t) — iJi(O)]^), the mean 
square change since the start in the measured specific energy 
of the particles. It can be seen that the value of this quantity 
rises roughly linearly with time, as HB90 found, indicating 
that the values are changing through a diffusive process. I 
fit a straight line to the last 90 values (i.e. ignoring the first 
11 , where the rise is often a little steeper), and henceforth 


report only this fitted slope as the mean square change of E 
per dynamical time to, and its associated uncertainty. 

The changes presented in Figure [3] are averages over all 
the particles. Naturally, the rms changes are larger for those 
particles whose orbit periods are shorter. As this trend is 
quite gradual in the inhomogeneous models and non-existent 
in the uniform models, the global averages I report in all 
cases are representative of those at intermediate radii where 
the bulk of the particles are located. 

The total angular momentum, L = |ij|, of each parti¬ 
cle is also conserved in a smooth potential, and its changes 
afford another measure of relaxation. I have found that the 
mean-squared changes in L also rise roughly linearly with 
time, but not quite as smoothly as the E changes in the left 
panel of Figure O Furthermore, mass segregation is still less 
evident in the angular momentum changes. I therefore focus 
on the E changes for the remainder of the paper. 


4.3 Mass segregation 

The right panel of Figure [3] gives the time evolution of the 
mean specific energy {Ei{t)) of the particles of the separate 
masses Q The rapid variations have the same sign for each 
particle species because they arise from variations in T+2W, 
as discussed in §3, that are related to the virial fluctuations 
(Fig. [21). They result from evolving potential changes seeded 
by shot-noise driven fluctuations among the particles. 

In addition to these fluctuations, the right panel of Fig¬ 
ure [3] displays a slow divergence of the mean energies of the 
heavy and light particles, which represents the gradual ex¬ 
change of energy that would in the long-run cause the heavy 
particles to settle to the centre and the lighter to populate 
the envelope of the model. The slopes of the straight lines 
fitted to these data should differ in magnitude by the ra¬ 
tio of the particle masses, since total energy conservation 
requires the energy lost by the heavy particles to be taken 


^ The mean energy measured from the simulation, (Ei{0)) ^ 
—0.30, is somewhat higher than that expected from the DF, 
(E) fa —0.356, because the potential well in the simulation deter¬ 
mined from the selected particles, which are only 74% of the total 
mass, is not as deep as the analytic potential of the untruncated 
model. 
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up by the light. In practice, this symmetry is imperfect be¬ 
cause of the large measurement uncertainties. Henceforth, I 
report only the slopes, multiplied by /r, and their statistical 
uncertainties. 

In systems of gravitating particles, heavy particles 
lose energy to light part icles through dynamical friction 
llBinnev fc Tremainell200§ l. Physically, a particle is braked 
by the attraction from its wake - i.e. the response of the 
surrounding sea of particles to its motion. In a system of 
equal mass particles, the frictional drag on each pa rticle is 
balan ced, on average, by the scattering accelerations (iHenonI 
Il973f l and no secular changes occur. But differences between 
the braking and acceleration terms cause mass segregation 
when particle masses are unequal. 

Were energy exchange between particles of the different 
mass species the only source of relaxation, the trends in the 
right panel of Figure [3] would not be strongly masked by 
short term oscillations. That short-term changes are larger 
than the gradual diverging trend is evidence that relaxation, 
measured in the left-hand panel, is being driven mostly by 
the collective oscillations of the model discussed above that 
are independent of particle masses. 

4.4 Measures of relaxation rate 

Figure |4] summarizes the relaxation measurements from all 
27 simulations. The three rows show the three different mass 
models; within each panel there are three different numbers 
of particles, and for each case, the evolution was computed 
by three separate force determination methods. The mea¬ 
sured rates from the heavy particles are shown in red, while 
those from the light particles are marked in green, and the 
uncertainties in the slopes are indicated by the error bars, 
that are often too short to be visible. The different force de¬ 
termination methods used are distinguished by the different 
symbol types as indicated in the figure caption, and offset 
horizontally from each other for clarity even though the par¬ 
ticle numbers are the same. These conventions are the same 
in the left and right panels. 

The left panels show the energy diffusion rate, dehned 
as 

{Emr^ (s) 

in units of . The adopted values for (i5i(0)), which are 
—0.3, —0.45, —0.9 for the Hernquist, Plummer, and uniform 
sphere respectively, are the same for all nine simulations 
with each model. The straight dotted line has slope —0.5 
while the slope of the dashed line is —1; these lines are not 
hts to the data and are for comparison only. The right panels 
show Nd{Ei{t))/dt, the factor N for each sub-population is 
included for clarity - in reality, the slopes roughly decrease 
as 1/N. 

There are many conclusions to be drawn from these 
data. 

First, the differences between the energy diffusion rates 
(left panels) for the heavy (red) and light (green) parti¬ 
cles within one simulation are generally small. This is a 
further indication of the dominance of coherent potential 
variations, which cause deflections that are independent of 
particle mass, in driving this measure of relaxation. 

Second, the energy diffusion rate declines roughly as 


in the uniform sphere, whereas in the inhomogeneous 
models it declines more or less with the expected N~^ de¬ 
pendence for collisional relaxation at fixed resolution (see 
§1). If two-body effects were dominant in all cases, the vari¬ 
ation with N should be the same in all three mass models. 

Third, the rates of energy exchange shown in the right 
hand panels generally have opposite signs, and vary roughly 
as since they are approximately constant when mul¬ 

tiplied by N. The values from the highest N experiments 
are quite uncertain, because the decreasing variation in the 
mean as N rises is masked by short-term changes (see the 
right panel of Fig. 0 - i.e. the trends decrease into the 
noise. This is particularly problematic for the uniform sphere 
(bottom right panel), where the potential changes associ¬ 
ated with pulsations of the model dominate over the energy 
exchange rate at the highest N, making an accurate mea¬ 
surement over the time interval simulated impossible. 

Fourth, the different methods used to compute the grav¬ 
itational field yield broadly similar behaviour in all three 
mass models, and the variation of the rates with N in both 
the left and right panels is similar for each method in each 
of the different models. 

Fifth, the energy diffusion rate (left panels) is consis¬ 
tently lower, but by a factor < 3, for the SFP method (di¬ 
amonds) than for the tree code (crosses) and the spherical 
grid (circles) in the Hernquist and Plummer models. This 
systematic difference i s also present in the variations of L. 
iHernauist fc Qstrikerl lll992fl also reported a slightly lower 
energy diffusion rate when they used a held method, al¬ 
though it was unclear whether the more rapid diffusion in 
their tree code, for example, resulted mainly from the mild 
disequilbrium of their initial model caused by gravity soft¬ 
ening. However, the mass segregation rate (right panels), 
which is a more direct consequence of two-body scattering, 
is no less for SFP than for the other methods. 

Sixth, relaxation in the uniform sphere probably is com¬ 
pletely dominated by the oscillatory modes excited by the 
shot noise in the particle distribution. The initial ampli¬ 
tudes of the modes should scale as and the decline 

in relaxation rate with approximately this dependence in the 
uniform sphere is a direct indication of their dominance in 
this case. Collective modes are present in all models, but 
the different A'^-dependence in the inhomogeneous models is 
probably because modes are rapidly damped in those cases. 

Seventh, short-range gravity softening appears to have 
little effect on the relaxation rate. There is no explicit soft¬ 
ening in SFP, and the only smoothing on the spherical grid 
is linear interpolation in radius. Force softening is explicit 
in the tree code and implicit in the C3D grid, where forces 
are slightly sharper CeH = 0.036a while e = 0.05a for the 
BHT code for the uniform sphere. Yet the relaxation rates 
scarcely differ between the various codes. This emphasizes 
the dominance of distant encounters, enhanced by collective 
oscillations, in driving relaxation and that softening’s only 
value is to avoid large accelerations during close encounters 
between particles, which would require short time steps to 
integrate the motions accurately. 

Once again, the quantities shown in Figure |4] are aver¬ 
ages over all the particles, and these conclusions therefore 
apply to the particles in the bulk of the model. The small 
number of more tightly bound particles have the largest en¬ 
ergy changes, but by factors of only a few, and therefore have 
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Figure 4. Summary of results from 27 simulations. The left column shows the A^-dependence of the diffusion rate defined in eq. lO and 
the right column shows Nd(Ei(t))/dt, with green points for the light particles and red for the nine times heavier particles. The top row 
is for the Hernquist model, the middle row for the Plummer model and the bottom row is for the uniform sphere. Crosses are from the 
tree code (BHT), circles from the spherical grid (S3D), diamonds from the field method (SFP) which was used for the Hernquist and 
Plummer models only, while squares are for a cubic Cartesian grid (C3D) which was used for the uniform sphere only. The error bars 
show dicr uncertainties in the slopes, and the symbols from the different methods have been slightly shifted horizontally for clarity even 
though the same total numbers of particles, N, were used. Note that the vertical range is expanded by a factor of 10 for the points at 
the highest N in the lower right panel. 


© 2015 RAS, MNRAS 000, [TH9] 





















8 J. A. Sellwood 


little overall effect on the global average. Note also that their 
rate of energy change should be reckoned on the time-scale 
of their shorter orbit periods. 


5 DISCUSSION 

5.1 Discrepancy with IWeinberg &: Kat j (I2007air5l 

The finding stated in point four above is in strong disagree¬ 
ment with that of lWeinberg &: Kat j ll2007alf^ . who argued 
that the relaxation time-scale is orders of magnitude longer 
in field methods. It should be noted that their analysis ad¬ 
dressed the more limited problem of resonant exchanges be¬ 
tween particles and a perturbing bar potential, not the more 
general relaxation studied her e. Howeve r , a di screpancy re¬ 
mains to be explained since ISellwoodI (l2008l l. in a study 
that addressed their specific prediction, also found that field 
methods were not superior to grid methods. 

Their calculation is highly technical but, in essence, 
they began by separating large- from small-scale noise and 
calculated separate contributions to the orbit deflections 
from both these “components” of the noise. Rather surpris¬ 
ingly, they found that deflections due to large-scale noise 
were much weaker than those from small-scale noise. They 
then argued that direct codes possess noise on both scales 
while field methods are affected only by large-scale noise, 
which led them to conclude that relaxation should be far 
slower in field methods. However, it is well-known that noise 
is present on all scales and, since every decade of impact 

J iarameter contributes equ ally to the Coulomb logarithm 
Binnev fc Tremain^ |200^ ~) , there can be no clear distinc¬ 
tion between large- and small-scale noise. This also implies 
that relaxation from large-scale noise should be comparable 
to that from small-scale noise. 

5.2 More on collective effects 

IVasilievI (120151 ') found good agreement between experimen¬ 
tal relaxation rates and theoretical predictions from Fokker- 
Planck diffusion coefficients, suggesting that self-gravity of 
the collective modes does not boost the relaxation rate sig¬ 
nificantly. Furthermore, his predicted rate for the Plummer 
sphere is in good agreement with the value I find (Fig. 21 
middle row, left panel). 

However, such studies usually consider inhomogeneous 
models. The anomalous results presented here for the uni¬ 
form sphere show that relaxation is greatly enhanced by 
collective modes in this un usual case where all orbits have 
the same period. Note that iRauch fc Tremaine! (Il996l l also 
found enhanced relaxation through collective effects in their 
study of a system of low-mass particles in near Keplerian 
motion about a central mass. In their case, it was important 
that particles with the same L had the same period. 

5.3 Mild dependence on method 

The reason for the slightly lower diffusion rate in the field 
method (point five above) is unclear. Both the SFP and 
S3D methods employ an expansion in surface harmonics up 
to Imax = 8 to capture any angular variations, and the prin¬ 
cipal difference between them is radial resolution: the SFP 


method employs 11 radial functions (0 ^ n Si 10), while the 
S3D uses a grid of 100 radial shells. 

It therefore seemed that the marginally higher energy 
diffusion rate when the S3D grid was used was because 
that grid could support more oscillatory modes. However, 
increasing the number of radial functions used in the SFP 
method to 51, in order to enable the field method also to 
support more collective modes, led to little change in either 
measure of the relaxation rate, confounding this theory! Ad¬ 
ditional experiments with /max = 4 and Zmax = 16 caused 
very slight changes in the measured relaxation rates in the 
expected sense, but the systematic discrepancy between the 
two methods persisted. 

While the puzzle remains, the negative outcome of these 
tests is further evidence of the unimportance of self-gravity 
for collective oscillations in inhomogeneous models. 


6 CONCLUSIONS 

The main result he r e con firms that previously found by 
iHernauist fc BarnesI (Il990l l and by iHernauist fc Ostrikerl 
I I 992 II that the rate of relaxation in N-body simulations 
that aspire to be collisionless is very largely independent of 
the method used to compute the gravitational field. This 
is especially true when the relaxation rate is assayed as the 
rate of energy exchange between particles of different masses 
- see the right hand panels of Fig. [d] This result is physi¬ 
cally reasonable, since relaxation is dominated by distant 
encounters and any method that correctly yields the field 
from distant particles must faithfully include their stochas¬ 
tic contributions. 

The rate of relaxation arises from at least two distin¬ 
guishable sources: a slow diffusion of the integrals caused 
by coherent potential oscillations of the system that are 
largely independent of particle mass, and mass segregation 
that is more directly caused by two-body encounters. The 
good agreement (§5.2) between Fokker-Planck diffusion and 
N-body experiment in the inhomogeneous models indicates 
that self-gravity of the collective modes adds little to the re¬ 
laxation rate. The rate of energy exchange between particles 
of differing masses becomes more difficult to measure as N 
rises because mean energies scatter with potential variations 
that scale as , while the trends in the mean energies 

of the different mass species diverge as N~^. 

The different N-dependence in the uniform sphere (bot¬ 
tom left panel of Fig. 0 is clear evidence that collective 
modes do dominate the relaxation in this special case. This 
model differs from the other two by having a harmonic po¬ 
tential throughout in which the orbit frequencies of all par¬ 
ticles are the same; this difference permits undamped col¬ 
lective oscillations, whereas collective modes are damped in 
the other cases. 

The slightly lower energy diffusion rates that result from 
use of the field method (diamonds in top left and middle 
left panels of Fig. 0 is bought at a high price. The leading 
term in the basis used for both cases was a perfect match 
to the equilibrium model, and the same basis would be less 
suited were the density to evolve, or were It used for any 
other model, requiring more terms to yield the correct to¬ 
tal potential. Thus field methods lack the versatility to fol¬ 
low arbitrary changes to the distribution of mass within a 
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model unless the expansion is taken to higher order, and 
their slightly better relaxation rate can be achieved in any 
of the more general methods simply by employing a few 
times mo re particles. A real advant age of field methods, fea¬ 
tured by IWdnberg &: Katd (l2007al la). is that they offer an 
elegant comparison of simulations with perturbation theory 
when computing first order changes to an equilibrium model 
that could be unstable or externally perturbed. 

The energy diffusion time-scale, the inverse of the rate 
defined in eq. (O and plotted in Fig. IH is > 10® dynam¬ 
ical times (i.e. 10^^ yr for the suggested scaling) for only 
= 4 X 10® particles in the inhomogeneous models (it is 
shorter for the exotic case of the uniform sphere). By this 
measure, relaxation times in simulations of inhomogeneous 
3D models using any valid code are already significantly 
longer than the ages of the galaxies being simulated, and 
relaxation considerations alone do not require much larger 
numbers of particles. 
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